Construct Single-Hierarchical P/NBD Model for Short Timeframe Synthetic Data

In this workbook we construct the non-hierarchical P/NBD models on the synthetic data with the longer timeframe.

1 Load and Construct Datasets

We start by modelling the P/NBD model using our synthetic datasets before we try to model real-life data.

1.1 Load Short Time-frame Synthetic Data

We now want to load the short time-frame synthetic data.

Code

customer_cohortdata_tbl <- read_rds("data/synthdata_shortframe_cohort_tbl.rds")
customer_cohortdata_tbl |> glimpse()

Rows: 50,000
Columns: 4
$ customer_id    <chr> "SFC202001_0001", "SFC202001_0002", "SFC202001_0003", "…
$ cohort_qtr     <chr> "2020 Q1", "2020 Q1", "2020 Q1", "2020 Q1", "2020 Q1", …
$ cohort_ym      <chr> "2020 01", "2020 01", "2020 01", "2020 01", "2020 01", …
$ first_tnx_date <date> 2020-01-01, 2020-01-01, 2020-01-01, 2020-01-01, 2020-0…

Code

customer_simparams_tbl  <- read_rds("data/synthdata_shortframe_simparams_tbl.rds")
customer_simparams_tbl |> glimpse()

Rows: 50,000
Columns: 9
$ customer_id     <chr> "SFC202001_0001", "SFC202001_0002", "SFC202001_0003", …
$ cohort_qtr      <chr> "2020 Q1", "2020 Q1", "2020 Q1", "2020 Q1", "2020 Q1",…
$ cohort_ym       <chr> "2020 01", "2020 01", "2020 01", "2020 01", "2020 01",…
$ first_tnx_date  <date> 2020-01-01, 2020-01-01, 2020-01-01, 2020-01-01, 2020-…
$ customer_lambda <dbl> 0.12213805, 0.29987747, 0.31504009, 0.03856001, 0.1881…
$ customer_mu     <dbl> 0.127118566, 0.096184402, 0.052334526, 0.204708842, 0.…
$ customer_tau    <dbl> 29.5956480, 10.9437199, 1.6938450, 2.6798108, 48.75206…
$ customer_amtmn  <dbl> 46.2371662, 111.0425353, 45.8870891, 43.0754249, 10.93…
$ customer_amtcv  <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …

Code

customer_transactions_tbl <- read_rds("data/synthdata_shortframe_transactions_tbl.rds")
customer_transactions_tbl |> glimpse()

Rows: 296,756
Columns: 4
$ customer_id   <chr> "SFC202001_0021", "SFC202001_0007", "SFC202001_0008", "S…
$ tnx_timestamp <dttm> 2020-01-01 00:34:40, 2020-01-01 02:41:45, 2020-01-01 04…
$ invoice_id    <chr> "T20200101-0001", "T20200101-0002", "T20200101-0003", "T…
$ tnx_amount    <dbl> 20.60, 1.83, 0.07, 30.03, 112.09, 45.85, 71.83, 166.33, …

1.2 Load Derived Datasets

Code

id_1000  <- read_rds("data/shortsynth_id_1000.rds")
id_5000  <- read_rds("data/shortsynth_id_5000.rds")
id_10000 <- read_rds("data/shortsynth_id_10000.rds")

fit_1000_data_tbl  <- read_rds("data/shortsynth_fit_1000_data_tbl.rds")
fit_10000_data_tbl <- read_rds("data/shortsynth_fit_10000_data_tbl.rds")

customer_fit_stats_tbl    <- fit_1000_data_tbl
customer_summarystats_tbl <- read_rds("data/shortsynth_customer_summarystats_tbl.rds")

obs_fitdata_tbl   <- read_rds("data/shortsynth_obs_fitdata_tbl.rds")
obs_validdata_tbl <- read_rds("data/shortsynth_obs_validdata_tbl.rds")

customer_subset_id <- id_1000

1.3 Load Subset Data

We also want to construct our data subsets for the purposes of speeding up our valuations.

Code

customer_fit_subset_tbl <- obs_fitdata_tbl |>
  filter(customer_id %in% customer_subset_id)

customer_fit_subset_tbl |> glimpse()

Rows: 1,000
Columns: 6
$ customer_id    <chr> "SFC202001_0016", "SFC202001_0027", "SFC202001_0094", "…
$ first_tnx_date <dttm> 2020-01-01 10:27:01, 2020-01-01 11:35:35, 2020-01-03 1…
$ last_tnx_date  <dttm> 2021-12-30 03:24:56, 2020-02-24 23:10:44, 2020-01-03 1…
$ tnx_count      <dbl> 80, 4, 0, 2, 7, 1, 6, 4, 2, 58, 0, 1, 15, 8, 37, 3, 0, …
$ t_x            <dbl> 104.1009825, 7.7832485, 0.0000000, 1.6433938, 47.262213…
$ T_cal          <dbl> 104.3664, 104.3596, 104.0277, 104.0590, 104.0572, 103.8…

Code

customer_valid_subset_tbl <- obs_validdata_tbl |>
  filter(customer_id %in% customer_subset_id)

customer_valid_subset_tbl |> glimpse()

Rows: 1,000
Columns: 3
$ customer_id       <chr> "SFC202001_0016", "SFC202001_0027", "SFC202001_0094"…
$ tnx_count         <dbl> 36, 0, 0, 0, 0, 0, 0, 0, 0, 26, 0, 0, 6, 0, 15, 0, 0…
$ tnx_last_interval <dbl> 51.45780, NA, NA, NA, NA, NA, NA, NA, NA, 50.27866, …

We now use these datasets to set the start and end dates for our various validation methods.

Code

dates_lst <- read_rds("data/shortsynth_simulation_dates.rds")

use_fit_start_date <- dates_lst$use_fit_start_date
use_fit_end_date   <- dates_lst$use_fit_end_date

use_valid_start_date <- dates_lst$use_valid_start_date
use_valid_end_date   <- dates_lst$use_valid_end_date

Finally, we need to set our directories where we save our Stan code and the model outputs.

Code

stan_modeldir <- "stan_models"
stan_codedir  <-   "stan_code"

2 Fit First Hierarchical Lambda Model

Our first hierarchical model puts a hierarchical prior around the mean of our population \(\lambda\) - lambda_mn.

Once again we use a Gamma prior for it.

2.1 Compile and Fit Stan Model

We now compile this model using CmdStanR.

Code

pnbd_onehierlambda_stanmodel <- cmdstan_model(
  "stan_code/pnbd_onehier_lambda.stan",
  include_paths =   stan_codedir,
  pedantic      =           TRUE,
  dir           =  stan_modeldir
  )

We then use this compiled model with our data to produce a fit of the data.

Code

stan_modelname <- "pnbd_init_short_onehierlambda1"
stanfit_prefix <- str_c("fit_", stan_modelname)
stanfit_seed   <- stanfit_seed + 1

stanfit_object_file <- glue("data/{stanfit_prefix}_stanfit.rds")


stan_data_lst <- customer_fit_stats_tbl |>
  select(customer_id, x, t_x, T_cal) |>
  compose_data(
    hier_lambda_mn_p1 = 0.25,
    hier_lambda_mn_p2 = 1,

    lambda_cv = 1.00,
    
    mu_mn     = 0.10,
    mu_cv     = 1.00,
    )

if(!file_exists(stanfit_object_file)) {
  pnbd_init_short_onehierlambda1_stanfit <- pnbd_onehierlambda_stanmodel$sample(
    data            =                stan_data_lst,
    chains          =                            4,
    iter_warmup     =                          500,
    iter_sampling   =                          500,
    seed            =                 stanfit_seed,
    save_warmup     =                         TRUE,
    output_dir      =                stan_modeldir,
    output_basename =               stanfit_prefix,
    )
  
  pnbd_init_short_onehierlambda1_stanfit$save_object(stanfit_object_file, compress = "gzip")

} else {
  pnbd_init_short_onehierlambda1_stanfit <- read_rds(stanfit_object_file)
}

Running MCMC with 4 chains, at most 12 in parallel...

Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 4 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 3 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 4 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 3 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 2 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 3 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 4 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 2 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 4 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 3 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 2 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 4 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 3 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 4 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 3 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 4 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 3 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 4 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 3 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 4 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 3 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 4 finished in 19.0 seconds.
Chain 1 finished in 19.1 seconds.
Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 3 finished in 19.4 seconds.
Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 2 finished in 19.8 seconds.

All 4 chains finished successfully.
Mean chain execution time: 19.3 seconds.
Total execution time: 20.0 seconds.

Code

pnbd_init_short_onehierlambda1_stanfit$summary()

# A tibble: 3,003 × 10
   variable       mean   median      sd     mad       q5      q95  rhat ess_bulk
   <chr>         <num>    <num>   <num>   <num>    <num>    <num> <num>    <num>
 1 lp__       -1.48e+4 -1.48e+4 36.0    36.2    -1.48e+4 -1.47e+4 1.01      581.
 2 lambda_mn   2.54e-1  2.54e-1  0.0109  0.0105  2.37e-1  2.73e-1 1.00     2226.
 3 lambda[1]   7.48e-1  7.45e-1  0.0827  0.0808  6.18e-1  8.88e-1 1.00     5877.
 4 lambda[2]   3.62e-1  3.39e-1  0.174   0.164   1.28e-1  6.76e-1 1.00     4006.
 5 lambda[3]   1.41e-1  8.02e-2  0.169   0.0919  5.43e-3  4.64e-1 0.999    2514.
 6 lambda[4]   4.12e-1  3.64e-1  0.270   0.233   8.53e-2  9.43e-1 1.00     3601.
 7 lambda[5]   1.40e-1  1.34e-1  0.0518  0.0497  6.86e-2  2.36e-1 1.00     3077.
 8 lambda[6]   1.39e-1  1.14e-1  0.105   0.0856  1.95e-2  3.41e-1 0.999    3156.
 9 lambda[7]   4.29e-1  4.08e-1  0.175   0.168   1.86e-1  7.59e-1 1.00     4181.
10 lambda[8]   3.14e-1  2.94e-1  0.144   0.134   1.18e-1  5.79e-1 1.00     3262.
# ℹ 2,993 more rows
# ℹ 1 more variable: ess_tail <num>

We have some basic HMC-based validity statistics we can check.

Code

pnbd_init_short_onehierlambda1_stanfit$cmdstan_diagnose()

Processing csv files: /home/rstudio/btydwork/stan_models/fit_pnbd_init_short_onehierlambda1-1.csvWarning: non-fatal error reading adaptation data
, /home/rstudio/btydwork/stan_models/fit_pnbd_init_short_onehierlambda1-2.csvWarning: non-fatal error reading adaptation data
, /home/rstudio/btydwork/stan_models/fit_pnbd_init_short_onehierlambda1-3.csvWarning: non-fatal error reading adaptation data
, /home/rstudio/btydwork/stan_models/fit_pnbd_init_short_onehierlambda1-4.csvWarning: non-fatal error reading adaptation data


Checking sampler transitions treedepth.
Treedepth satisfactory for all transitions.

Checking sampler transitions for divergences.
No divergent transitions found.

Checking E-BFMI - sampler transitions HMC potential energy.
E-BFMI satisfactory.

Effective sample size satisfactory.

Split R-hat values satisfactory all parameters.

Processing complete, no problems detected.

2.2 Visual Diagnostics of the Sample Validity

Now that we have a sample from the posterior distribution we need to create a few different visualisations of the diagnostics.

Code

parameter_subset <- c(
  "lambda_mn",
  "lambda[1]", "lambda[2]", "lambda[3]", "lambda[4]",
  "mu[1]",     "mu[2]",     "mu[3]",     "mu[4]"
  )

pnbd_init_short_onehierlambda1_stanfit$draws(inc_warmup = FALSE) |>
  mcmc_trace(pars = parameter_subset) +
  expand_limits(y = 0) +
  labs(
    x = "Iteration",
    y = "Value",
    title = "Traceplot of Sample of Lambda and Mu Values"
    ) +
  theme(axis.text.x = element_text(size = 10))

We also check \(N_{eff}\) as a quick diagnostic of the fit.

Code

pnbd_init_short_onehierlambda1_stanfit |>
  neff_ratio(pars = c("lambda", "mu")) |>
  as.numeric() |>
  mcmc_neff() +
    ggtitle("Plot of Parameter Effective Sample Sizes")

2.3 Assess the Model

As we intend to run the same logic to assess each of our models, we have combined all this logic into a single function run_model_assessment, to run the simulations and combine the datasets.

Code

pnbd_stanfit <- pnbd_init_short_onehierlambda1_stanfit |>
  recover_types(customer_fit_stats_tbl)

pnbd_init_short_onehierlambda1_assess_data_lst <- run_model_assessment(
  model_stanfit    = pnbd_stanfit,
  insample_tbl     = customer_fit_subset_tbl,
  fit_label        = "pnbd_init_short_onehierlambda1",
  fit_end_dttm     = use_fit_end_date     |> as.POSIXct(),
  valid_start_dttm = use_valid_start_date |> as.POSIXct(),
  valid_end_dttm   = use_valid_end_date   |> as.POSIXct(),
  sim_seed         = 4210
  )

pnbd_init_short_onehierlambda1_assess_data_lst |> glimpse()

List of 5
 $ model_fit_index_filepath     : 'glue' chr "data/pnbd_init_short_onehierlambda1_assess_fit_index_tbl.rds"
 $ model_valid_index_filepath   : 'glue' chr "data/pnbd_init_short_onehierlambda1_assess_valid_index_tbl.rds"
 $ model_simstats_filepath      : 'glue' chr "data/pnbd_init_short_onehierlambda1_assess_model_simstats_tbl.rds"
 $ model_fit_simstats_filepath  : 'glue' chr "data/pnbd_init_short_onehierlambda1_assess_fit_simstats_tbl.rds"
 $ model_valid_simstats_filepath: 'glue' chr "data/pnbd_init_short_onehierlambda1_assess_valid_simstats_tbl.rds"

2.3.1 Check In-Sample Data Validation

We first check the model against the in-sample data.

Code

simdata_tbl <- pnbd_init_short_onehierlambda1_assess_data_lst |>
  use_series(model_fit_simstats_filepath) |>
  read_rds()

insample_plots_lst <- create_model_assessment_plots(
  obsdata_tbl = obs_fitdata_tbl,
  simdata_tbl = simdata_tbl
  )

insample_plots_lst$multi_plot |> print()

Code

insample_plots_lst$total_plot |> print()

Code

insample_plots_lst$quant_plot |> print()

This fit looks reasonable and appears to capture most of the aspects of the data used to fit it. Given that this is a synthetic dataset, this is not surprising, but at least we appreciate that our model is valid.

2.3.2 Check Out-of-Sample Data Validation

We now repeat for the out-of-sample data.

Code

simdata_tbl <- pnbd_init_short_onehierlambda1_assess_data_lst |>
  use_series(model_valid_simstats_filepath) |>
  read_rds()

outsample_plots_lst <- create_model_assessment_plots(
  obsdata_tbl = obs_validdata_tbl,
  simdata_tbl = simdata_tbl
  )

outsample_plots_lst$multi_plot |> print()

Code

outsample_plots_lst$total_plot |> print()

Code

outsample_plots_lst$quant_plot |> print()

As for our short time frame data, overall our model is working well.

3 Fit Second Hierarchical Lambda Model

In this model, we are going with a broadly similar model but we are instead using a different mean for our hierarchical prior.

3.1 Fit Stan Model

We now want to fit the model to our data using this alternative prior for lambda_mn.

Code

stan_modelname <- "pnbd_init_short_onehierlambda2"
stanfit_prefix <- str_c("fit_", stan_modelname)
stanfit_seed   <- stanfit_seed + 1

stanfit_object_file <- glue("data/{stanfit_prefix}_stanfit.rds")


stan_data_lst <- customer_fit_stats_tbl |>
  select(customer_id, x, t_x, T_cal) |>
  compose_data(
    hier_lambda_mn_p1 = 0.50,
    hier_lambda_mn_p2 = 1,

    lambda_cv = 1.00,
    
    mu_mn     = 0.10,
    mu_cv     = 1.00,
    )

if(!file_exists(stanfit_object_file)) {
  pnbd_init_short_onehierlambda2_stanfit <- pnbd_onehierlambda_stanmodel$sample(
    data            =                stan_data_lst,
    chains          =                            4,
    iter_warmup     =                          500,
    iter_sampling   =                          500,
    seed            =                 stanfit_seed,
    save_warmup     =                         TRUE,
    output_dir      =                stan_modeldir,
    output_basename =               stanfit_prefix,
    )
  
  pnbd_init_short_onehierlambda2_stanfit$save_object(stanfit_object_file, compress = "gzip")

} else {
  pnbd_init_short_onehierlambda2_stanfit <- read_rds(stanfit_object_file)
}

Running MCMC with 4 chains, at most 12 in parallel...

Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 3 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 4 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 2 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 3 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 4 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 2 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 3 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 4 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 2 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 4 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 3 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 4 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 3 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 3 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 4 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 3 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 4 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 3 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 4 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 3 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 4 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 1 finished in 20.9 seconds.
Chain 2 finished in 21.0 seconds.
Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 3 finished in 21.5 seconds.
Chain 4 finished in 21.5 seconds.

All 4 chains finished successfully.
Mean chain execution time: 21.2 seconds.
Total execution time: 21.7 seconds.

Code

pnbd_init_short_onehierlambda2_stanfit$summary()

# A tibble: 3,003 × 10
   variable       mean   median      sd     mad       q5      q95  rhat ess_bulk
   <chr>         <num>    <num>   <num>   <num>    <num>    <num> <num>    <num>
 1 lp__       -1.48e+4 -1.48e+4 36.7    36.0    -1.48e+4 -1.47e+4 1.01      555.
 2 lambda_mn   2.55e-1  2.55e-1  0.0106  0.0106  2.37e-1  2.72e-1 1.00     2724.
 3 lambda[1]   7.48e-1  7.44e-1  0.0840  0.0892  6.17e-1  8.89e-1 1.00     4375.
 4 lambda[2]   3.60e-1  3.32e-1  0.174   0.163   1.29e-1  6.78e-1 1.00     3579.
 5 lambda[3]   1.43e-1  8.40e-2  0.168   0.0966  5.72e-3  4.75e-1 1.00     2990.
 6 lambda[4]   4.12e-1  3.50e-1  0.268   0.241   9.84e-2  9.10e-1 1.00     2879.
 7 lambda[5]   1.41e-1  1.34e-1  0.0522  0.0516  6.66e-2  2.36e-1 0.999    4473.
 8 lambda[6]   1.43e-1  1.12e-1  0.114   0.0893  2.13e-2  3.60e-1 1.00     3956.
 9 lambda[7]   4.28e-1  4.03e-1  0.174   0.165   1.84e-1  7.53e-1 1.00     4873.
10 lambda[8]   3.17e-1  2.98e-1  0.150   0.139   1.20e-1  5.91e-1 1.00     3602.
# ℹ 2,993 more rows
# ℹ 1 more variable: ess_tail <num>

We have some basic HMC-based validity statistics we can check.

Code

pnbd_init_short_onehierlambda2_stanfit$cmdstan_diagnose()

Processing csv files: /home/rstudio/btydwork/stan_models/fit_pnbd_init_short_onehierlambda2-1.csvWarning: non-fatal error reading adaptation data
, /home/rstudio/btydwork/stan_models/fit_pnbd_init_short_onehierlambda2-2.csvWarning: non-fatal error reading adaptation data
, /home/rstudio/btydwork/stan_models/fit_pnbd_init_short_onehierlambda2-3.csvWarning: non-fatal error reading adaptation data
, /home/rstudio/btydwork/stan_models/fit_pnbd_init_short_onehierlambda2-4.csvWarning: non-fatal error reading adaptation data


Checking sampler transitions treedepth.
Treedepth satisfactory for all transitions.

Checking sampler transitions for divergences.
No divergent transitions found.

Checking E-BFMI - sampler transitions HMC potential energy.
E-BFMI satisfactory.

Effective sample size satisfactory.

Split R-hat values satisfactory all parameters.

Processing complete, no problems detected.

3.2 Visual Diagnostics of the Sample Validity

Now that we have a sample from the posterior distribution we need to create a few different visualisations of the diagnostics.

Code

parameter_subset <- c(
  "lambda_mn",
  "lambda[1]", "lambda[2]", "lambda[3]", "lambda[4]",
  "mu[1]",     "mu[2]",     "mu[3]",     "mu[4]"
  )

pnbd_init_short_onehierlambda2_stanfit$draws(inc_warmup = FALSE) |>
  mcmc_trace(pars = parameter_subset) +
  expand_limits(y = 0) +
  labs(
    x = "Iteration",
    y = "Value",
    title = "Traceplot of Sample of Lambda and Mu Values"
    ) +
  theme(axis.text.x = element_text(size = 10))

We also check \(N_{eff}\) as a quick diagnostic of the fit.

Code

pnbd_init_short_onehierlambda2_stanfit |>
  neff_ratio(pars = c("lambda", "mu")) |>
  as.numeric() |>
  mcmc_neff() +
    ggtitle("Plot of Parameter Effective Sample Sizes")

3.3 Assess the Model

As we intend to run the same logic to assess each of our models, we have combined all this logic into a single function run_model_assessment, to run the simulations and combine the datasets.

Code

pnbd_stanfit <- pnbd_init_short_onehierlambda2_stanfit |>
  recover_types(customer_fit_stats_tbl)

pnbd_init_short_onehierlambda2_assess_data_lst <- run_model_assessment(
  model_stanfit    = pnbd_stanfit,
  insample_tbl     = customer_fit_subset_tbl,
  fit_label        = "pnbd_init_short_onehierlambda2",
  fit_end_dttm     = use_fit_end_date     |> as.POSIXct(),
  valid_start_dttm = use_valid_start_date |> as.POSIXct(),
  valid_end_dttm   = use_valid_end_date   |> as.POSIXct(),
  sim_seed         = 4210
  )

pnbd_init_short_onehierlambda2_assess_data_lst |> glimpse()

List of 5
 $ model_fit_index_filepath     : 'glue' chr "data/pnbd_init_short_onehierlambda2_assess_fit_index_tbl.rds"
 $ model_valid_index_filepath   : 'glue' chr "data/pnbd_init_short_onehierlambda2_assess_valid_index_tbl.rds"
 $ model_simstats_filepath      : 'glue' chr "data/pnbd_init_short_onehierlambda2_assess_model_simstats_tbl.rds"
 $ model_fit_simstats_filepath  : 'glue' chr "data/pnbd_init_short_onehierlambda2_assess_fit_simstats_tbl.rds"
 $ model_valid_simstats_filepath: 'glue' chr "data/pnbd_init_short_onehierlambda2_assess_valid_simstats_tbl.rds"

3.3.1 Check In-Sample Data Validation

We first check the model against the in-sample data.

Code

simdata_tbl <- pnbd_init_short_onehierlambda2_assess_data_lst |>
  use_series(model_fit_simstats_filepath) |>
  read_rds()

insample_plots_lst <- create_model_assessment_plots(
  obsdata_tbl = obs_fitdata_tbl,
  simdata_tbl = simdata_tbl
  )

insample_plots_lst$multi_plot |> print()

Code

insample_plots_lst$total_plot |> print()

Code

insample_plots_lst$quant_plot |> print()

This fit looks reasonable and appears to capture most of the aspects of the data used to fit it. Given that this is a synthetic dataset, this is not surprising, but at least we appreciate that our model is valid.

3.3.2 Check Out-of-Sample Data Validation

We now repeat for the out-of-sample data.

Code

simdata_tbl <- pnbd_init_short_onehierlambda1_assess_data_lst |>
  use_series(model_valid_simstats_filepath) |>
  read_rds()

outsample_plots_lst <- create_model_assessment_plots(
  obsdata_tbl = obs_validdata_tbl,
  simdata_tbl = simdata_tbl
  )

outsample_plots_lst$multi_plot |> print()

Code

outsample_plots_lst$total_plot |> print()

Code

outsample_plots_lst$quant_plot |> print()

As for our short time frame data, overall our model is working well.

4 Fit First Hierarchical Mu Model

We now construct the same hierarchical model but based around mu_mn.

4.1 Compile and Fit Stan Model

We compile this model using CmdStanR.

Code

pnbd_onehiermu_stanmodel <- cmdstan_model(
  "stan_code/pnbd_onehier_mu.stan",
  include_paths =   stan_codedir,
  pedantic      =           TRUE,
  dir           =  stan_modeldir
  )

We then use this compiled model with our data to produce a fit of the data.

Code

stan_modelname <- "pnbd_init_short_onehiermu1"
stanfit_prefix <- str_c("fit_", stan_modelname)
stanfit_seed   <- stanfit_seed + 1

stanfit_object_file <- glue("data/{stanfit_prefix}_stanfit.rds")


stan_data_lst <- customer_fit_stats_tbl |>
  select(customer_id, x, t_x, T_cal) |>
  compose_data(
    hier_mu_mn_p1 = 0.50,
    hier_mu_mn_p2 = 1.00,

    lambda_mn     = 0.25,
    lambda_cv     = 1.00,

    mu_cv         = 1.00

    )

if(!file_exists(stanfit_object_file)) {
  pnbd_init_short_onehiermu1_stanfit <- pnbd_onehiermu_stanmodel$sample(
    data            =                stan_data_lst,
    chains          =                            4,
    iter_warmup     =                          500,
    iter_sampling   =                          500,
    seed            =                 stanfit_seed,
    save_warmup     =                         TRUE,
    output_dir      =                stan_modeldir,
    output_basename =               stanfit_prefix,
    )
  
  pnbd_init_short_onehiermu1_stanfit$save_object(stanfit_object_file, compress = "gzip")

} else {
  pnbd_init_short_onehiermu1_stanfit <- read_rds(stanfit_object_file)
}

Running MCMC with 4 chains, at most 12 in parallel...

Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 3 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 4 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 3 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 2 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 4 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 2 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 3 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 4 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 3 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 2 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 4 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 3 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 4 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 3 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 4 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 3 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 4 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 3 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 4 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 3 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 4 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 3 finished in 22.4 seconds.
Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 2 finished in 22.6 seconds.
Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 1 finished in 22.8 seconds.
Chain 4 finished in 22.7 seconds.

All 4 chains finished successfully.
Mean chain execution time: 22.6 seconds.
Total execution time: 22.9 seconds.

Code

pnbd_init_short_onehiermu1_stanfit$summary()

# A tibble: 3,003 × 10
   variable       mean   median      sd     mad       q5      q95  rhat ess_bulk
   <chr>         <num>    <num>   <num>   <num>    <num>    <num> <num>    <num>
 1 lp__       -1.39e+4 -1.39e+4 3.56e+1 3.61e+1 -1.39e+4 -1.38e+4  1.00     638.
 2 mu_mn       1.04e-1  1.04e-1 7.51e-3 7.59e-3  9.25e-2  1.17e-1  1.00    1027.
 3 lambda[1]   7.46e-1  7.44e-1 8.12e-2 8.30e-2  6.19e-1  8.83e-1  1.01    5017.
 4 lambda[2]   3.63e-1  3.33e-1 1.79e-1 1.69e-1  1.33e-1  7.06e-1  1.00    4549.
 5 lambda[3]   1.44e-1  8.45e-2 1.67e-1 9.52e-2  4.96e-3  4.75e-1  1.00    2211.
 6 lambda[4]   4.09e-1  3.50e-1 2.64e-1 2.21e-1  9.73e-2  9.11e-1  1.00    3636.
 7 lambda[5]   1.39e-1  1.32e-1 5.09e-2 4.95e-2  6.82e-2  2.31e-1  1.00    3497.
 8 lambda[6]   1.42e-1  1.11e-1 1.14e-1 9.21e-2  1.89e-2  3.68e-1  1.00    2946.
 9 lambda[7]   4.28e-1  4.05e-1 1.73e-1 1.62e-1  1.87e-1  7.52e-1  1.00    3712.
10 lambda[8]   3.19e-1  2.90e-1 1.54e-1 1.43e-1  1.18e-1  6.11e-1  1.01    5044.
# ℹ 2,993 more rows
# ℹ 1 more variable: ess_tail <num>

We have some basic HMC-based validity statistics we can check.

Code

pnbd_init_short_onehiermu1_stanfit$cmdstan_diagnose()

Processing csv files: /home/rstudio/btydwork/stan_models/fit_pnbd_init_short_onehiermu1-1.csvWarning: non-fatal error reading adaptation data
, /home/rstudio/btydwork/stan_models/fit_pnbd_init_short_onehiermu1-2.csvWarning: non-fatal error reading adaptation data
, /home/rstudio/btydwork/stan_models/fit_pnbd_init_short_onehiermu1-3.csvWarning: non-fatal error reading adaptation data
, /home/rstudio/btydwork/stan_models/fit_pnbd_init_short_onehiermu1-4.csvWarning: non-fatal error reading adaptation data


Checking sampler transitions treedepth.
Treedepth satisfactory for all transitions.

Checking sampler transitions for divergences.
No divergent transitions found.

Checking E-BFMI - sampler transitions HMC potential energy.
E-BFMI satisfactory.

Effective sample size satisfactory.

Split R-hat values satisfactory all parameters.

Processing complete, no problems detected.

4.2 Visual Diagnostics of the Sample Validity

Now that we have a sample from the posterior distribution we need to create a few different visualisations of the diagnostics.

Code

parameter_subset <- c(
  "mu_mn",
  "lambda[1]", "lambda[2]", "lambda[3]", "lambda[4]",
  "mu[1]",     "mu[2]",     "mu[3]",     "mu[4]"
  )

pnbd_init_short_onehiermu1_stanfit$draws(inc_warmup = FALSE) |>
  mcmc_trace(pars = parameter_subset) +
  expand_limits(y = 0) +
  labs(
    x = "Iteration",
    y = "Value",
    title = "Traceplot of Sample of Lambda and Mu Values"
    ) +
  theme(axis.text.x = element_text(size = 10))

We also check \(N_{eff}\) as a quick diagnostic of the fit.

Code

pnbd_init_short_onehiermu1_stanfit |>
  neff_ratio(pars = c("lambda", "mu")) |>
  as.numeric() |>
  mcmc_neff() +
    ggtitle("Plot of Parameter Effective Sample Sizes")

4.3 Assess the Model

As we intend to run the same logic to assess each of our models, we have combined all this logic into a single function run_model_assessment, to run the simulations and combine the datasets.

Code

pnbd_stanfit <- pnbd_init_short_onehiermu1_stanfit |>
  recover_types(customer_fit_stats_tbl)

pnbd_init_short_onehiermu1_assess_data_lst <- run_model_assessment(
  model_stanfit    = pnbd_stanfit,
  insample_tbl     = customer_fit_subset_tbl,
  fit_label        = "pnbd_init_short_onehiermu1",
  fit_end_dttm     = use_fit_end_date     |> as.POSIXct(),
  valid_start_dttm = use_valid_start_date |> as.POSIXct(),
  valid_end_dttm   = use_valid_end_date   |> as.POSIXct(),
  sim_seed         = 4210
  )

pnbd_init_short_onehiermu1_assess_data_lst |> glimpse()

List of 5
 $ model_fit_index_filepath     : 'glue' chr "data/pnbd_init_short_onehiermu1_assess_fit_index_tbl.rds"
 $ model_valid_index_filepath   : 'glue' chr "data/pnbd_init_short_onehiermu1_assess_valid_index_tbl.rds"
 $ model_simstats_filepath      : 'glue' chr "data/pnbd_init_short_onehiermu1_assess_model_simstats_tbl.rds"
 $ model_fit_simstats_filepath  : 'glue' chr "data/pnbd_init_short_onehiermu1_assess_fit_simstats_tbl.rds"
 $ model_valid_simstats_filepath: 'glue' chr "data/pnbd_init_short_onehiermu1_assess_valid_simstats_tbl.rds"

4.3.1 Check In-Sample Data Validation

We first check the model against the in-sample data.

Code

simdata_tbl <- pnbd_init_short_onehiermu1_assess_data_lst |>
  use_series(model_fit_simstats_filepath) |>
  read_rds()

insample_plots_lst <- create_model_assessment_plots(
  obsdata_tbl = obs_fitdata_tbl,
  simdata_tbl = simdata_tbl
  )

insample_plots_lst$multi_plot |> print()

Code

insample_plots_lst$total_plot |> print()

Code

insample_plots_lst$quant_plot |> print()

This fit looks reasonable and appears to capture most of the aspects of the data used to fit it. Given that this is a synthetic dataset, this is not surprising, but at least we appreciate that our model is valid.

4.3.2 Check Out-of-Sample Data Validation

We now repeat for the out-of-sample data.

Code

simdata_tbl <- pnbd_init_short_onehierlambda1_assess_data_lst |>
  use_series(model_valid_simstats_filepath) |>
  read_rds()

outsample_plots_lst <- create_model_assessment_plots(
  obsdata_tbl = obs_validdata_tbl,
  simdata_tbl = simdata_tbl
  )

outsample_plots_lst$multi_plot |> print()

Code

outsample_plots_lst$total_plot |> print()

Code

outsample_plots_lst$quant_plot |> print()

As for our short time frame data, overall our model is working well.

5 Fit Second Hierarchical Mu Model

In this model, we are going with a broadly similar model but we are instead using a different mean for our hierarchical prior.

5.1 Fit Stan Model

We now want to fit the model to our data using this alternative prior for lambda_mn.

Code

stan_modelname <- "pnbd_init_short_onehiermu2"
stanfit_prefix <- str_c("fit_", stan_modelname)
stanfit_seed   <- stanfit_seed + 1

stanfit_object_file <- glue("data/{stanfit_prefix}_stanfit.rds")


stan_data_lst <- customer_fit_stats_tbl |>
  select(customer_id, x, t_x, T_cal) |>
  compose_data(
    hier_mu_mn_p1 = 0.25,
    hier_mu_mn_p2 = 1.00,

    lambda_mn     = 0.25,
    lambda_cv     = 1.00,

    mu_cv         = 1.00
    )

if(!file_exists(stanfit_object_file)) {
  pnbd_init_short_onehiermu2_stanfit <- pnbd_onehiermu_stanmodel$sample(
    data            =                stan_data_lst,
    chains          =                            4,
    iter_warmup     =                          500,
    iter_sampling   =                          500,
    seed            =                 stanfit_seed,
    save_warmup     =                         TRUE,
    output_dir      =                stan_modeldir,
    output_basename =               stanfit_prefix,
    )
  
  pnbd_init_short_onehiermu2_stanfit$save_object(stanfit_object_file, compress = "gzip")

} else {
  pnbd_init_short_onehiermu2_stanfit <- read_rds(stanfit_object_file)
}

Running MCMC with 4 chains, at most 12 in parallel...

Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 3 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 4 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 2 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 3 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 4 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 2 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 3 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 4 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 2 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 3 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 4 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 3 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 4 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 3 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 4 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 3 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 4 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 3 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 4 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 3 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 4 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 1 finished in 17.8 seconds.
Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 2 finished in 18.0 seconds.
Chain 3 finished in 17.9 seconds.
Chain 4 finished in 17.9 seconds.

All 4 chains finished successfully.
Mean chain execution time: 17.9 seconds.
Total execution time: 18.1 seconds.

Code

pnbd_init_short_onehiermu2_stanfit$summary()

# A tibble: 3,003 × 10
   variable       mean   median      sd     mad       q5      q95  rhat ess_bulk
   <chr>         <num>    <num>   <num>   <num>    <num>    <num> <num>    <num>
 1 lp__       -1.39e+4 -1.39e+4 3.58e+1 3.59e+1 -1.39e+4 -1.38e+4  1.00     611.
 2 mu_mn       1.04e-1  1.04e-1 7.25e-3 7.33e-3  9.27e-2  1.17e-1  1.00     923.
 3 lambda[1]   7.50e-1  7.46e-1 8.43e-2 8.73e-2  6.19e-1  8.95e-1  1.01    6602.
 4 lambda[2]   3.59e-1  3.32e-1 1.71e-1 1.60e-1  1.34e-1  6.73e-1  1.00    3151.
 5 lambda[3]   1.40e-1  8.54e-2 1.66e-1 9.32e-2  4.62e-3  4.72e-1  1.00    2595.
 6 lambda[4]   4.10e-1  3.48e-1 2.64e-1 2.28e-1  9.23e-2  9.24e-1  1.00    3038.
 7 lambda[5]   1.41e-1  1.34e-1 5.04e-2 5.01e-2  7.00e-2  2.30e-1  1.00    4122.
 8 lambda[6]   1.40e-1  1.10e-1 1.08e-1 8.42e-2  2.25e-2  3.70e-1  1.00    3596.
 9 lambda[7]   4.33e-1  4.06e-1 1.77e-1 1.60e-1  1.91e-1  7.51e-1  1.00    4319.
10 lambda[8]   3.16e-1  2.89e-1 1.51e-1 1.46e-1  1.17e-1  6.03e-1  1.00    3420.
# ℹ 2,993 more rows
# ℹ 1 more variable: ess_tail <num>

We have some basic HMC-based validity statistics we can check.

Code

pnbd_init_short_onehiermu2_stanfit$cmdstan_diagnose()

Processing csv files: /home/rstudio/btydwork/stan_models/fit_pnbd_init_short_onehiermu2-1.csvWarning: non-fatal error reading adaptation data
, /home/rstudio/btydwork/stan_models/fit_pnbd_init_short_onehiermu2-2.csvWarning: non-fatal error reading adaptation data
, /home/rstudio/btydwork/stan_models/fit_pnbd_init_short_onehiermu2-3.csvWarning: non-fatal error reading adaptation data
, /home/rstudio/btydwork/stan_models/fit_pnbd_init_short_onehiermu2-4.csvWarning: non-fatal error reading adaptation data


Checking sampler transitions treedepth.
Treedepth satisfactory for all transitions.

Checking sampler transitions for divergences.
No divergent transitions found.

Checking E-BFMI - sampler transitions HMC potential energy.
E-BFMI satisfactory.

Effective sample size satisfactory.

Split R-hat values satisfactory all parameters.

Processing complete, no problems detected.

5.2 Visual Diagnostics of the Sample Validity

Now that we have a sample from the posterior distribution we need to create a few different visualisations of the diagnostics.

Code

parameter_subset <- c(
  "mu_mn",
  "lambda[1]", "lambda[2]", "lambda[3]", "lambda[4]",
  "mu[1]",     "mu[2]",     "mu[3]",     "mu[4]"
  )

pnbd_init_short_onehiermu2_stanfit$draws(inc_warmup = FALSE) |>
  mcmc_trace(pars = parameter_subset) +
  expand_limits(y = 0) +
  labs(
    x = "Iteration",
    y = "Value",
    title = "Traceplot of Sample of Lambda and Mu Values"
    ) +
  theme(axis.text.x = element_text(size = 10))

We also check \(N_{eff}\) as a quick diagnostic of the fit.

Code

pnbd_init_short_onehiermu2_stanfit |>
  neff_ratio(pars = c("lambda", "mu")) |>
  as.numeric() |>
  mcmc_neff() +
    ggtitle("Plot of Parameter Effective Sample Sizes")

5.3 Assess the Model

As we intend to run the same logic to assess each of our models, we have combined all this logic into a single function run_model_assessment, to run the simulations and combine the datasets.

Code

pnbd_stanfit <- pnbd_init_short_onehiermu2_stanfit |>
  recover_types(customer_fit_stats_tbl)

pnbd_init_short_onehiermu2_assess_data_lst <- run_model_assessment(
  model_stanfit    = pnbd_stanfit,
  insample_tbl     = customer_fit_subset_tbl,
  fit_label        = "pnbd_init_short_onehiermu2",
  fit_end_dttm     = use_fit_end_date     |> as.POSIXct(),
  valid_start_dttm = use_valid_start_date |> as.POSIXct(),
  valid_end_dttm   = use_valid_end_date   |> as.POSIXct(),
  sim_seed         = 4210
  )

pnbd_init_short_onehiermu2_assess_data_lst |> glimpse()

List of 5
 $ model_fit_index_filepath     : 'glue' chr "data/pnbd_init_short_onehiermu2_assess_fit_index_tbl.rds"
 $ model_valid_index_filepath   : 'glue' chr "data/pnbd_init_short_onehiermu2_assess_valid_index_tbl.rds"
 $ model_simstats_filepath      : 'glue' chr "data/pnbd_init_short_onehiermu2_assess_model_simstats_tbl.rds"
 $ model_fit_simstats_filepath  : 'glue' chr "data/pnbd_init_short_onehiermu2_assess_fit_simstats_tbl.rds"
 $ model_valid_simstats_filepath: 'glue' chr "data/pnbd_init_short_onehiermu2_assess_valid_simstats_tbl.rds"

5.3.1 Check In-Sample Data Validation

We first check the model against the in-sample data.

Code

simdata_tbl <- pnbd_init_short_onehiermu2_assess_data_lst |>
  use_series(model_fit_simstats_filepath) |>
  read_rds()

insample_plots_lst <- create_model_assessment_plots(
  obsdata_tbl = obs_fitdata_tbl,
  simdata_tbl = simdata_tbl
  )

insample_plots_lst$multi_plot |> print()

Code

insample_plots_lst$total_plot |> print()

Code

insample_plots_lst$quant_plot |> print()

This fit looks reasonable and appears to capture most of the aspects of the data used to fit it. Given that this is a synthetic dataset, this is not surprising, but at least we appreciate that our model is valid.

5.3.2 Check Out-of-Sample Data Validation

We now repeat for the out-of-sample data.

Code

simdata_tbl <- pnbd_init_short_onehiermu1_assess_data_lst |>
  use_series(model_valid_simstats_filepath) |>
  read_rds()

outsample_plots_lst <- create_model_assessment_plots(
  obsdata_tbl = obs_validdata_tbl,
  simdata_tbl = simdata_tbl
  )

outsample_plots_lst$multi_plot |> print()

Code

outsample_plots_lst$total_plot |> print()

Code

outsample_plots_lst$quant_plot |> print()

As for our short time frame data, overall our model is working well.

6 Compare Model Outputs

We have looked at each of the models individually, but it is also worth looking at each of the models as a group.

Code

calculate_simulation_statistics <- function(file_rds) {
  simdata_tbl <- read_rds(file_rds)
  
  multicount_cust_tbl <- simdata_tbl |>
    filter(sim_tnx_count > 0) |>
    count(draw_id, name = "multicust_count")
  
  totaltnx_data_tbl <- simdata_tbl |>
    count(draw_id, wt = sim_tnx_count, name = "simtnx_count")
  
  simstats_tbl <- multicount_cust_tbl |>
    inner_join(totaltnx_data_tbl, by = "draw_id")
  
  return(simstats_tbl)
}

Code

obs_fit_customer_count <- customer_fit_subset_tbl |>
  filter(tnx_count > 0) |>
  nrow()

obs_valid_customer_count <- customer_valid_subset_tbl |>
  filter(tnx_count > 0) |>
  nrow()

obs_fit_total_count <- customer_fit_subset_tbl |>
  pull(tnx_count) |>
  sum()

obs_valid_total_count <- customer_valid_subset_tbl |>
  pull(tnx_count) |>
  sum()


obs_stats_tbl <- tribble(
  ~assess_type, ~name,               ~obs_value,
  "fit",        "multicust_count",   obs_fit_customer_count,
  "fit",        "simtnx_count",      obs_fit_total_count,
  "valid",      "multicust_count",   obs_valid_customer_count,
  "valid",      "simtnx_count",      obs_valid_total_count
  )


model_assess_tbl <- dir_ls("data", regexp = "pnbd_init_short_(one|fixed).*_assess_.*simstats") |>
  enframe(name = NULL, value = "file_path") |>
  filter(str_detect(file_path, "_assess_model_", negate = TRUE)) |>
  mutate(
    model_label = str_replace(file_path, "data/pnbd_init_short_(.*?)_assess_.*", "\\1"),
    assess_type = if_else(str_detect(file_path, "_assess_fit_"), "fit", "valid"),
    
    sim_data = map(
      file_path, calculate_simulation_statistics,
      
      .progress = "calculate_simulation_statistics"
      )
    )

model_assess_tbl |> glimpse()

Rows: 14
Columns: 4
$ file_path   <fs::path> "data/pnbd_init_short_fixed1_assess_fit_simstats_tbl.…
$ model_label <chr> "fixed1", "fixed1", "fixed2", "fixed2", "fixed3", "fixed3"…
$ assess_type <chr> "fit", "valid", "fit", "valid", "fit", "valid", "fit", "va…
$ sim_data    <list> [<tbl_df[2000 x 3]>], [<tbl_df[2000 x 3]>], [<tbl_df[2000…

Code

model_assess_summstat_tbl <- model_assess_tbl |>
  select(model_label, assess_type, sim_data) |>
  unnest(sim_data) |>
  pivot_longer(
    cols = !c(model_label, assess_type, draw_id)
    ) |>
  group_by(model_label, assess_type, name) |>
  summarise(
    .groups = "drop",
    
    mean_val = mean(value),
    p10 = quantile(value, 0.10),
    p25 = quantile(value, 0.25),
    p50 = quantile(value, 0.50),
    p75 = quantile(value, 0.75),
    p90 = quantile(value, 0.90)
    )

model_assess_summstat_tbl |> glimpse()

Rows: 28
Columns: 9
$ model_label <chr> "fixed1", "fixed1", "fixed1", "fixed1", "fixed2", "fixed2"…
$ assess_type <chr> "fit", "fit", "valid", "valid", "fit", "fit", "valid", "va…
$ name        <chr> "multicust_count", "simtnx_count", "multicust_count", "sim…
$ mean_val    <dbl> 619.8015, 4314.2750, 178.1055, 1822.4220, 681.0555, 3200.6…
$ p10         <dbl> 601.0, 4045.9, 168.0, 1655.0, 664.0, 3010.9, 127.0, 822.9,…
$ p25         <dbl> 611.00, 4168.00, 173.00, 1737.00, 672.00, 3093.00, 132.00,…
$ p50         <dbl> 620.0, 4310.0, 178.0, 1820.5, 681.0, 3196.0, 137.0, 931.0,…
$ p75         <dbl> 629, 4457, 183, 1908, 690, 3301, 141, 993, 540, 5183, 540,…
$ p90         <dbl> 637.0, 4593.1, 188.0, 1984.1, 698.0, 3397.0, 146.0, 1046.0…

Code

#! echo: TRUE

ggplot(model_assess_summstat_tbl) +
  geom_errorbar(
    aes(x = model_label, ymin = p10, ymax = p90), width = 0
    ) +
  geom_errorbar(
    aes(x = model_label, ymin = p25, ymax = p75), width = 0, linewidth = 3
    ) +
  geom_hline(
    aes(yintercept = obs_value),
    data = obs_stats_tbl, colour = "red"
    ) +
  scale_y_continuous(labels = label_comma()) +
  expand_limits(y = 0) +
  facet_wrap(
    vars(assess_type, name), scale = "free_y"
    ) +
  labs(
    x = "Model",
    y = "Count",
    title = "Comparison Plot for the Different Models"
    ) +
  theme(
    axis.text.x = element_text(angle = 20, vjust = 0.5, size = 8)
    )

6.1 Write Assessment Data to Disk

We now want to save the assessment data to disk.

Code

model_assess_tbl |> write_rds("data/assess_data_pnbd_init_short_onehier_tbl.rds")

7 R Environment

Code

options(width = 120L)
sessioninfo::session_info()

─ Session info ───────────────────────────────────────────────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.3.1 (2023-06-16)
 os       Ubuntu 22.04.3 LTS
 system   x86_64, linux-gnu
 ui       X11
 language (EN)
 collate  en_US.UTF-8
 ctype    en_US.UTF-8
 tz       Europe/Dublin
 date     2023-10-19
 pandoc   3.1.1 @ /usr/local/bin/ (via rmarkdown)

─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────
 package        * version    date (UTC) lib source
 abind            1.4-5      2016-07-21 [1] RSPM (R 4.3.1)
 arrayhelpers     1.1-0      2020-02-04 [1] RSPM (R 4.3.1)
 backports        1.4.1      2021-12-13 [1] RSPM (R 4.3.0)
 base64enc        0.1-3      2015-07-28 [1] RSPM (R 4.3.0)
 bayesplot      * 1.10.0     2022-11-16 [1] RSPM (R 4.3.1)
 bit              4.0.5      2022-11-15 [1] RSPM (R 4.3.0)
 bit64            4.0.5      2020-08-30 [1] RSPM (R 4.3.0)
 bridgesampling   1.1-2      2021-04-16 [1] RSPM (R 4.3.1)
 brms           * 2.20.4     2023-09-25 [1] RSPM (R 4.3.1)
 Brobdingnag      1.2-9      2022-10-19 [1] RSPM (R 4.3.1)
 cachem           1.0.8      2023-05-01 [1] RSPM (R 4.3.0)
 callr            3.7.3      2022-11-02 [1] RSPM (R 4.3.0)
 checkmate        2.2.0      2023-04-27 [1] RSPM (R 4.3.1)
 cli              3.6.1      2023-03-23 [1] RSPM (R 4.3.0)
 cmdstanr       * 0.6.0.9000 2023-10-17 [1] Github (stan-dev/cmdstanr@a13c798)
 coda             0.19-4     2020-09-30 [1] RSPM (R 4.3.1)
 codetools        0.2-19     2023-02-01 [2] CRAN (R 4.3.1)
 colorspace       2.1-0      2023-01-23 [1] RSPM (R 4.3.0)
 colourpicker     1.3.0      2023-08-21 [1] RSPM (R 4.3.1)
 conflicted     * 1.2.0      2023-02-01 [1] RSPM (R 4.3.1)
 cowplot        * 1.1.1      2020-12-30 [1] RSPM (R 4.3.1)
 crayon           1.5.2      2022-09-29 [1] RSPM (R 4.3.0)
 crosstalk        1.2.0      2021-11-04 [1] RSPM (R 4.3.1)
 data.table       1.14.8     2023-02-17 [1] RSPM (R 4.3.0)
 digest           0.6.33     2023-07-07 [1] RSPM (R 4.3.0)
 directlabels   * 2023.8.25  2023-09-01 [1] RSPM (R 4.3.1)
 distributional   0.3.2      2023-03-22 [1] RSPM (R 4.3.1)
 dplyr          * 1.1.3      2023-09-03 [1] RSPM (R 4.3.0)
 DT               0.30       2023-10-05 [1] RSPM (R 4.3.1)
 dygraphs         1.1.1.6    2018-07-11 [1] RSPM (R 4.3.1)
 ellipsis         0.3.2      2021-04-29 [1] RSPM (R 4.3.0)
 evaluate         0.21       2023-05-05 [1] RSPM (R 4.3.0)
 fansi            1.0.4      2023-01-22 [1] RSPM (R 4.3.0)
 farver           2.1.1      2022-07-06 [1] RSPM (R 4.3.0)
 fastmap          1.1.1      2023-02-24 [1] RSPM (R 4.3.0)
 forcats        * 1.0.0      2023-01-29 [1] RSPM (R 4.3.0)
 fs             * 1.6.3      2023-07-20 [1] RSPM (R 4.3.1)
 furrr          * 0.3.1      2022-08-15 [1] RSPM (R 4.3.1)
 future         * 1.33.0     2023-07-01 [1] RSPM (R 4.3.1)
 generics         0.1.3      2022-07-05 [1] RSPM (R 4.3.0)
 ggdist           3.3.0      2023-05-13 [1] RSPM (R 4.3.1)
 ggplot2        * 3.4.3      2023-08-14 [1] RSPM (R 4.3.0)
 globals          0.16.2     2022-11-21 [1] RSPM (R 4.3.1)
 glue           * 1.6.2      2022-02-24 [1] RSPM (R 4.3.0)
 gridExtra        2.3        2017-09-09 [1] RSPM (R 4.3.1)
 gtable           0.3.4      2023-08-21 [1] RSPM (R 4.3.0)
 gtools           3.9.4      2022-11-27 [1] RSPM (R 4.3.1)
 hms              1.1.3      2023-03-21 [1] RSPM (R 4.3.0)
 htmltools        0.5.6      2023-08-10 [1] RSPM (R 4.3.0)
 htmlwidgets      1.6.2      2023-03-17 [1] RSPM (R 4.3.0)
 httpuv           1.6.11     2023-05-11 [1] RSPM (R 4.3.0)
 igraph           1.5.1      2023-08-10 [1] RSPM (R 4.3.1)
 inline           0.3.19     2021-05-31 [1] RSPM (R 4.3.1)
 jsonlite         1.8.7      2023-06-29 [1] RSPM (R 4.3.0)
 knitr            1.44       2023-09-11 [1] RSPM (R 4.3.0)
 labeling         0.4.3      2023-08-29 [1] RSPM (R 4.3.0)
 later            1.3.1      2023-05-02 [1] RSPM (R 4.3.0)
 lattice          0.21-8     2023-04-05 [2] CRAN (R 4.3.1)
 lifecycle        1.0.3      2022-10-07 [1] RSPM (R 4.3.0)
 listenv          0.9.0      2022-12-16 [1] RSPM (R 4.3.1)
 loo              2.6.0      2023-03-31 [1] RSPM (R 4.3.1)
 lubridate      * 1.9.2      2023-02-10 [1] RSPM (R 4.3.0)
 magrittr       * 2.0.3      2022-03-30 [1] RSPM (R 4.3.0)
 markdown         1.9        2023-09-30 [1] RSPM (R 4.3.1)
 Matrix           1.5-4.1    2023-05-18 [2] CRAN (R 4.3.1)
 matrixStats      1.0.0      2023-06-02 [1] RSPM (R 4.3.1)
 memoise          2.0.1      2021-11-26 [1] RSPM (R 4.3.0)
 mime             0.12       2021-09-28 [1] RSPM (R 4.3.0)
 miniUI           0.1.1.1    2018-05-18 [1] RSPM (R 4.3.0)
 munsell          0.5.0      2018-06-12 [1] RSPM (R 4.3.0)
 mvtnorm          1.2-3      2023-08-25 [1] RSPM (R 4.3.1)
 nlme             3.1-162    2023-01-31 [2] CRAN (R 4.3.1)
 parallelly       1.36.0     2023-05-26 [1] RSPM (R 4.3.1)
 pillar           1.9.0      2023-03-22 [1] RSPM (R 4.3.0)
 pkgbuild         1.4.2      2023-06-26 [1] RSPM (R 4.3.0)
 pkgconfig        2.0.3      2019-09-22 [1] RSPM (R 4.3.0)
 plyr             1.8.9      2023-10-02 [1] RSPM (R 4.3.1)
 posterior      * 1.4.1      2023-03-14 [1] RSPM (R 4.3.1)
 prettyunits      1.2.0      2023-09-24 [1] RSPM (R 4.3.0)
 processx         3.8.2      2023-06-30 [1] RSPM (R 4.3.0)
 promises         1.2.1      2023-08-10 [1] RSPM (R 4.3.0)
 ps               1.7.5      2023-04-18 [1] RSPM (R 4.3.0)
 purrr          * 1.0.2      2023-08-10 [1] RSPM (R 4.3.0)
 quadprog         1.5-8      2019-11-20 [1] RSPM (R 4.3.1)
 QuickJSR         1.0.6      2023-09-12 [1] RSPM (R 4.3.1)
 R6               2.5.1      2021-08-19 [1] RSPM (R 4.3.0)
 Rcpp           * 1.0.11     2023-07-06 [1] RSPM (R 4.3.0)
 RcppParallel     5.1.7      2023-02-27 [1] RSPM (R 4.3.1)
 readr          * 2.1.4      2023-02-10 [1] RSPM (R 4.3.0)
 reshape2         1.4.4      2020-04-09 [1] RSPM (R 4.3.1)
 rlang          * 1.1.1      2023-04-28 [1] RSPM (R 4.3.0)
 rmarkdown        2.25       2023-09-18 [1] RSPM (R 4.3.0)
 rstan            2.26.23    2023-09-08 [1] RSPM (R 4.3.1)
 rstantools       2.3.1.1    2023-07-18 [1] RSPM (R 4.3.1)
 rstudioapi       0.15.0     2023-07-07 [1] RSPM (R 4.3.0)
 rsyslog        * 1.0.3      2023-05-08 [1] RSPM (R 4.3.1)
 scales         * 1.2.1      2022-08-20 [1] RSPM (R 4.3.0)
 sessioninfo      1.2.2      2021-12-06 [1] RSPM (R 4.3.1)
 shiny            1.7.5      2023-08-12 [1] RSPM (R 4.3.0)
 shinyjs          2.1.0      2021-12-23 [1] RSPM (R 4.3.1)
 shinystan        2.6.0      2022-03-03 [1] RSPM (R 4.3.1)
 shinythemes      1.2.0      2021-01-25 [1] RSPM (R 4.3.1)
 StanHeaders      2.26.28    2023-09-07 [1] RSPM (R 4.3.1)
 stringi          1.7.12     2023-01-11 [1] RSPM (R 4.3.0)
 stringr        * 1.5.0      2022-12-02 [1] RSPM (R 4.3.0)
 svUnit           1.0.6      2021-04-19 [1] RSPM (R 4.3.1)
 tensorA          0.36.2     2020-11-19 [1] RSPM (R 4.3.1)
 threejs          0.3.3      2020-01-21 [1] RSPM (R 4.3.1)
 tibble         * 3.2.1      2023-03-20 [1] RSPM (R 4.3.0)
 tidybayes      * 3.0.6      2023-08-12 [1] RSPM (R 4.3.1)
 tidyr          * 1.3.0      2023-01-24 [1] RSPM (R 4.3.0)
 tidyselect       1.2.0      2022-10-10 [1] RSPM (R 4.3.0)
 tidyverse      * 2.0.0      2023-02-22 [1] RSPM (R 4.3.0)
 timechange       0.2.0      2023-01-11 [1] RSPM (R 4.3.0)
 tzdb             0.4.0      2023-05-12 [1] RSPM (R 4.3.0)
 utf8             1.2.3      2023-01-31 [1] RSPM (R 4.3.0)
 vctrs            0.6.3      2023-06-14 [1] RSPM (R 4.3.0)
 vroom            1.6.3      2023-04-28 [1] RSPM (R 4.3.0)
 withr            2.5.0      2022-03-03 [1] RSPM (R 4.3.0)
 xfun             0.40       2023-08-09 [1] RSPM (R 4.3.0)
 xtable           1.8-4      2019-04-21 [1] RSPM (R 4.3.0)
 xts              0.13.1     2023-04-16 [1] RSPM (R 4.3.1)
 yaml             2.3.7      2023-01-23 [1] RSPM (R 4.3.0)
 zoo              1.8-12     2023-04-13 [1] RSPM (R 4.3.1)

 [1] /usr/local/lib/R/site-library
 [2] /usr/local/lib/R/library

──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────

Code

options(width = 80L)